#include "complexMat.h"
#include <cassert>

namespace moremath{

	Matd Matc::re(){
		Matd realMat(mNumRows, mNumCols, vl_0);
		for(int i=0; i<mNumRows; i++) for(int j=0; j<mNumCols; j++) realMat[i][j] = mMat[i][j].re();
		return realMat;
	}

	Matd Matc::im(){
		Matd imagMat(mNumRows, mNumCols, vl_0);
		for(int i=0; i<mNumRows; i++) for(int j=0; j<mNumCols; j++) imagMat[i][j] = mMat[i][j].im();
		return imagMat;

	}
	void Matc::setReal(Matd &_real){
		assert(_real.Rows()==mNumRows);
		assert(_real.Cols()==mNumCols);
		for(int i=0; i<mNumRows; i++) for(int j=0; j<mNumCols; j++) mMat[i][j].re() = _real[i][j];
	}
	void Matc::setImag(Matd &_imag){
		assert(_imag.Rows()==mNumRows);
		assert(_imag.Cols()==mNumCols);
		for(int i=0; i<mNumRows; i++) for(int j=0; j<mNumCols; j++) mMat[i][j].im() = _imag[i][j];
	}

	Vecc Matc::getCol(int _j){
		Vecc v(mNumRows);
		for(int i=0; i<mNumRows; i++) v[i] = mMat[i][_j];
		return v;
	}
	void Matc::setCol(int _j, Vecc _c){
		assert(_c.Elts()==mNumRows);
		for(int i=0; i<mNumRows; i++) mMat[i][_j] = _c[i];
	}

	void Matc::operator = (Matc _m){
		mNumRows = _m.Rows();
		mNumCols = _m.Cols();
		mMat.resize(mNumRows);
		for(int i=0; i<mNumRows; i++) mMat[i]= _m[i];
	}
	Matc Matc::operator + (Matc& _m){
		assert(_m.Rows()==mNumRows);
		assert(_m.Cols()==mNumCols);
		Matc mnew(mNumRows, mNumCols);
		for(int i=0; i<mNumRows; i++) mnew[i] = mMat[i] + _m[i];
		return mnew;
	}
	void Matc::operator += (Matc& _m){
		assert(_m.Rows()==mNumRows);
		assert(_m.Cols()==mNumCols);
		for(int i=0; i<mNumRows; i++) mMat[i] += _m[i];
	}
	Matc Matc::operator - (){
		Matc mnew(mNumRows, mNumCols);
		for(int i=0; i<mNumRows; i++) mnew[i] = -mMat[i];
		return mnew;
	}
	Matc Matc::operator * (double _s){
		Matc m(mNumRows, mNumCols);
		for(int i=0; i<mNumRows; i++) m[i]= mMat[i]*_s;
		return m;
	}
	Matc Matc::operator *(complex _q){
		Matc m(mNumRows, mNumCols);
		for(int i=0; i<mNumRows; i++) m[i]= mMat[i]*_q;
		return m;
	}

	void Matc::operator *= (double _s){
		for(int i=0; i<mNumRows; i++) mMat[i] *= _s;
	}

	Matc Matc::operator / (double _s){
		Matc m(mNumRows, mNumCols);
		for(int i=0; i<mNumRows; i++) m[i]= mMat[i]/_s;
		return m;
	}

	Vecc Matc::operator * (Vecc &_v){
		assert(mNumCols==_v.Elts());
		Vecc v(mNumRows);
		for(int i=0; i<mNumRows; i++) v[i] = cdot(mMat[i], _v);
		return v;
	}

	Vecc Matc::operator * (Vecd &_v){
		assert(mNumCols==_v.Elts());
		Vecc v(mNumRows);
		for(int i=0; i<mNumRows; i++) v[i] = cdot(mMat[i], _v);
		return v;
	}
	//Matc Matc::operator * (Matd &_m){
	//	assert(mNumCols==_m.Rows());		
	//	Matc mnew(mNumRows, _m.Cols());
	//	for(int i=0; i<mNumRows; i++){
	//		for(int j=0; j<_m.Cols(); j++){
	//			// dot the ith row of mat 1 with jth col of mat 2
	//			for(int k=0; k<mNumCols; k++){
	//				mnew[i][j]+=mMat[i][k]*_m[k][j];
	//			}
	//		}
	//	}
	//	return mnew;
	//}
	Matc Matc::operator * (Matc &_m){
		assert(mNumCols==_m.Rows());		
		Matc mnew(mNumRows, _m.Cols());
		for(int i=0; i<mNumRows; i++){
			for(int j=0; j<_m.Cols(); j++){
				// dot the ith row of mat 1 with jth col of mat 2
				for(int k=0; k<mNumCols; k++){
					mnew[i][j]+=mMat[i][k]*_m[k][j];
				}
			}
		}
		return mnew;
	}

	Matc Matc::conj(){
		Matc m(mNumRows, mNumCols);
		for(int i=0; i<mNumRows; i++) m[i]= mMat[i].conj();
		return m;
	}

}

moremath::Matc operator * (moremath::Matc &_m1, Matd &_m2){
	assert(_m1.Cols()==_m2.Rows());		
	moremath::Matc mnew(_m1.Rows(), _m2.Cols());
	for(int i=0; i<_m1.Rows(); i++){
		for(int j=0; j<_m2.Cols(); j++){
			// dot the ith row of mat 1 with jth col of mat 2
			for(int k=0; k<_m1.Cols(); k++){
				mnew[i][j]+=_m1[i][k]*_m2[k][j];
			}
		}
	}
	return mnew;
}

moremath::Matc operator * (Matd &_m1, moremath::Matc &_m2){
	assert(_m1.Cols()==_m2.Rows());		
	moremath::Matc mnew(_m1.Rows(), _m2.Cols());
	for(int i=0; i<_m1.Rows(); i++){
		for(int j=0; j<_m2.Cols(); j++){
			// dot the ith row of mat 1 with jth col of mat 2
			for(int k=0; k<_m1.Cols(); k++){
				mnew[i][j]+=_m2[k][j]*_m1[i][k];
			}
		}
	}
	return mnew;
}

moremath::Matc trans(moremath::Matc &_m){
	moremath::Matc mt(_m.Cols(), _m.Rows());
	for(int i=0; i<mt.Rows(); i++){
		mt[i] = _m.getCol(i);
	}
	return mt;
}

moremath::Vecc diag(moremath::Matc &_m){
	assert(_m.Rows()==_m.Cols());
	moremath::Vecc d(_m.Rows());
	for(int i=0; i<d.Elts(); i++) d[i] = _m[i][i];
	return d;
}
